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Abstract 

The  inclusion  of  anisotropic  surface  free  energy  and  anisotropic  linear  interface  ki- 
netics in  phase-field  models  is  studied  for  the  solidification  of  a pure  material.  The 
formulation  is  described  for  a two-dimensional  system  with  a smooth  crystal-melt  in- 
terface and  for  a surface  free  energy  that  varies  smoothly  with  orientation,  in  which 
case  a quite  general  dependence  of  the  surface  free  energy  and  kinetic  coefficient  on 
orientation  can  be  treated;  it  is  assumed  that  the  anisotropy  is  mUd  enough  that  miss- 
ing orientations  do  not  occur.  The  method  of  matched  asymptotic  expansions  is  used 
to  recover  the  appropriate  anisotropic  form  of  the  Gibbs-Thomson  equation  in  the 
sharp-interface  limit  in  which  the  width  of  the  diffuse  interface  is  thin  compared  to  its 
local  radius  of  curvature.  It  is  found  that  the  surface  free  energy  and  the  thickness  of 
the  diffuse  interface  have  the  same  anisotropy,  whereas  the  kinetic  coefficient  has  an 
anisotropy  characterized  by  the  product  of  the  interface  thickness  with  the  intrinsic 
mobility  of  the  phase  field. 
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1 Introduction 


Phase-field  models  [1,  2,  3,  4,  5,  6]  provide  a convenient  basis  for  the  numerical  solution  of 
complicated  solidification  problems.  In  a phase-field  model,  in  addition  to  the  customary 
energy  and/or  concentration  variables,  an  additional  variable,  the  phase  field,  if,  is  intro- 
duced to  label  exphcitly  the  liquid  and  sohd  phases.  The  phase  field  takes  on  a constant 
value  in  each  bulk  phase,  e.g.  ip  = 0 in  the  sohd  phase  and  c/?  — 1 in  the  liquid  phase. 
The  transformation  from  solid  to  liquid  occurs  over  a thin  transition  region  where  ip  varies 
smoothly  from  zero  to  one.  The  usual  thermodynamic  functions  describing  the  system  can 
then  be  modified  to  incorporate  gradient  energy  terms;  in  particular,  terms  proportional  to 
|V(,c>p  can  contribute  to  the  surface  excess  quantities  that  play  a fundamental  role  in  Gibbs’ 
formulation  of  surface  thermodynamics  [7].  In  this  sense,  phase- field  models  are  natural 
outgrowths  of  diffuse-interface  models  dating  back  to  work  by  Van  der  Waals  [8],  by  Cahn 
and  Allen  [9,  10],  and  by  Cahn  and  Hilhard  [11,  12].  From  a computational  viewpoint, 
phase- field  models  are  similar  in  some  ways  to  the  enthalpy  method  [13],  in  that  explicit 
tracking  of  the  solid-hquid  interface  is  avoided.  Phase-field  models,  however,  are  more  ver- 
satile than  enthalpy  methods  since  such  effects  as  undercoohng  of  the  melt  and  departures 
from  thermodynamic  equilibrium  at  the  interface  are  included  automatically  (see,  e.g.,  [14]). 

An  important  example  of  the  utility  of  phase-field  models  is  given  by  the  numerical  studies 
of  dendritic  growth  by  Kobayashi  [15,  16,  17,  18]  and  by  Wheeler,  Murray,  and  Schaefer  [19]. 
The  phase-field  treatments  seem  to  capture  successfully  a broad  variety  of  dendritic  growth 
phenomena,  including  the  correct  relation  between  Peclet  number  and  undercooling,  the 
emission  of  sidearms,  and  the  coarsening  behavior  of  sidearms  that  are  further  removed  from 
the  tip.  The  anisotropy  of  surface  free  energy  or  of  interface  kinetics  is  generally  thought  to 
play  a fundamental  role  in  the  dynamics  of  dendritic  growth  [20,  21]. 

Since  in  a phase-field  formulation  the  interface  is  diffuse,  the  proper  incorporation  of  sur- 
face free  energy  anisotropy  requires  careful  consideration.  Phase  field  models  with  anisotropy 
have  been  considered  previously  for  specific  choices  of  gradient  energy.  Caginalp  and  Fife 
[22,  23]  have  considered  models  in  which  the  isotropic  “square  gradient”  expression  is  replaced 
by  a more  general  quadratic  form  with  different  coefficients  in  each  coordinate  direction.  For 
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an  isothermal  system  this  leads  to  an  elliptical  equilibrium  shape.  In  order  to  obtain  more 
complicated  anisotropies,  Langer  [2]  proposed  the  addition  to  the  gradient  energy  of  terms 
involving  the  squares  of  higher  derivatives  of  the  phase  field,  and  gave  an  example  leading  to 
cubic  anisotropy.  Cahn  and  Kikuchi  [24]  have  considered  discrete  forms  of  diffuse  interface 
models,  and  have  also  included  anisotropic  effects  through  the  choice  of  nearest-neighbor  in- 
teractions (see  also  [25,  26]).  Both  Kobayashi  [17]  and  Wheeler  et  al.  [19]  include  anisotropy 
by  allowing  the  coefficient  of  the  gradient  energy  to  depend  on  the  local  orientation  of  the 
gradient  of  the  phase  field.  Early  numerical  calculations  were  performed  by  Smith  [27]  and 
by  Umantsev  et  al.  [28]  in  which  no  explicit  anisotropy  was  included  in  the  models;  rather, 
anisotropy  was  provided  imphcitly  by  the  underlying  grid  used  in  the  numerical  calculations. 

In  this  paper  we  present  an  asymptotic  analysis  in  the  sharp-interface  limit  of  the  model 
studied  by  Kobayashi  [17]  and  Wheeler  et  al.  [19],  including  an  anisotropic  mobility.  We 
consider  a two-dimensional  phase-field  description  for  the  solidification  of  a single  component 
material  of  uniform  density,  and  assume  that  there  is  no  convection  in  the  melt.  We  also 
assume  that  the  anisotropy  is  mild  enough  that  the  resulting  interface  shape  is  smooth.  In 
such  an  asymptotic  analysis,  the  width  of  the  transition  region  is  thin  compared  to  the  radius 
of  curvature  of  the  Hnes  = constant.  Similar  asymptotic  analyses  have  been  performed  for 
the  isotropic  case  [29,  30,  31,  32]  to  recover  the  boundary  condition 

- = Tm-  ^-tK  - T,  (1) 

fi  Lv 

that  relates  the  normal  velocity  of  the  interface  Vn  (considered  to  be  positive  for  the  formation 
of  sohd),  the  mean  curvature  of  the  interface,  K , and  the  interface  temperature,  T/;  here  fi  is 
the  interfacial  kinetic  coefficient,  Ta/  is  the  bulk  melting  point,  7 is  the  (isotropic)  interfacial 
free  energy,  and  Ly  is  the  latent  heat  of  fusion  per  unit  volume.  Our  goal  is  to  derive  from 
an  appropriate  phase-field  model  the  corresponding  anisotropic  form  of  this  equation  [33], 
which  in  two  dimensions  is 


(2) 


where  7 = 7(0)  denotes  the  dependence  of  the  surface  free  energy  on  the  local  interface 
orientation,  as  measured  by  the  angle  6 between  the  the  interface  normal  and  a given  crys- 
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tallographic  direction,  and  fi  = ii{9)  is  an  anisotropic  kinetic  coefficient  that  also  depends 
on  9',  here  'yee  denotes  the  second  derivative  of  the  function  'y{9). 


2 Isotropic  Formulation 


Of  prime  consideration  will  be  the  surface  free  energy,  7,  which  is  defined  in  terms  of  the 
surface  excess  of  the  Helmholtz  free  energy  of  the  system.  In  a phase-field  model,  the 
bulk  Helmoltz  free  energy  density,  /,  includes  dependence  on  the  phase  field,  9?,  so  that 
f = /(T,  (p),  where  T is  the  temperature.  The  free  energy  density  of  the  solid  is  then  f[T,  0) 
and  that  of  the  liquid  is  /(T,  1),  and  at  the  bulk  melting  point,  Tm,  the  two  are  equal. 
The  Helmholtz  free  energy  functional  that  gives  the  free  energy  of  an  isothermal  two-phase 
system  of  volume  V is  assumed  to  have  the  form  [34] 

^ = l{fiT,'p)+^f\V<p\^'^dV,  (3) 


where  the  constant  l-tiat  appears  in  the  gradient  energy  coefficient  has  units  of  energy  per 
unit  length.  If  we  adjust  the  bulk  free  energy  so  that  /(TmjO)  = /(3m,  1)  = 0,  then  for  a 
one- dimensional  system  with  (p  = (p{x)  the  surface  excess  free  energy  per  unit  area  is  given 

by  ^ 

1 = J |/(^M,  p)  + dx.  (4) 

For  an  isothermal  system,  an  evolution  equation  for  the  phase  field  is  often  postulated 
by  requiring  that  p evolve  so  as  to  minimize  that  is,  by  setting 


(5) 


where  tq  is  an  empirical  relaxation  coefficient,  whose  inverse  is  an  intrinsic  interfacial  mobil- 
ity. 


For  the  non-isothermal  case,  it  is  more  appropriate  to  start  with  an  entropy  functional 
for  the  system  [34,  35], 

S = (6) 

in  which  the  entropy  density  s depends  on  the  internal  energy  density  and  the  phase.  The 
parameter  ( that  appears  in  the  entropy  functional  may  be  related  to  the  parameter  ^0  of  tbe 
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Helmholtz  free  energy  functional  by  [34]-  Evolution  equations  for  the  temperature 

and  the  phase  field  may  be  derived  by  requiring  that  e and  ip  evolve  so  as  to  maximize  S\ 
for  example,  phase-field  equations  of  the  form 


de 

dt 


-V- 


) 


(7) 

(8) 


can  be  derived  in  this  manner  [34].  Here  Mt  is  proportional  to  the  thermal  conductivity, 
and  the  internal  energy  density  has  the  form 


e(T,  v)  = es(T)^  p{v)L{T)  = ez,(T)  + \p{v>)  - llL(r),  (9) 


where  L{Tm)  is  the  latent  heat  of  fusion  per  unit  volume  and  p{ip)  is  a smooth  function 
with  p(0)  = 0 and  p(l)  = 1,  so  that  e(T,  1)  — e(T,  0)  = — es{T)  — L{T).  The  energy 

equation  may  then  be  rewritten  in  the  form 


c{T,v>)^  + L{T)p'{^)^  = V ■ [kVT] , 


(10) 


where  we  have  introduced 


ciT,^)  = [1  - ^ + P(¥>)^, 


(11) 


which  is  an  interpolation  of  the  bulk  heat  capacities  per  unit  volume,  and  the  thermal 
conductivity,  k = the  thermal  conductivity  may  also  be  allowed  to  depend  on 

temperature  and  phase  in  a general  formulation.  The  function  Q{T)  is  given  by  the  expression 

[341 


Q{T)  = l 


T i(C)  _ L{Tm) 


I’m  C 


d(  = 


Th 


[T  - Tm]  + o(|r  - TmY) 


and  the  function  G{tp)  is  taken  to  be  a double-well  potential  of  the  form 

GM  = - <P)^ 


(12) 


(13) 


where  a is  a constant  that  determines  the  height  of  the  intermediate  maximum  in  the  double- 
well potential;  1/a  has  dimensions  of  energy  density  per  degree. 
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In  the  special  case  that  5e(T,0)/5T  = 5e(r,  1)/5T  = cq  and  L{T)  = Lq  are  constant, 
and  p{(p)  = ip,  Eqs.  (7)-(8)  become 


+ ffe 


(14) 


(15) 


which  is  essentially  equivalent  to  the  form  used  by  Langer  [1,  2]  and  Caginalp  [4].  A dis- 
advantage of  the  choice  p{<p)  = ip  is  that  the  roots  of  the  expression  Q[T)p'{ip)  — G'{ip) 
appearing  in  Eq.  (8)  that  determine  the  values  of  ip  in  the  bulk  phases  then  depend  on 
temperature;  choosing  a more  general  form  for  p{ip)  that  satisfies  p^(0)  = p'(l)  = 0 allows 
the  roots  ip  = Q and  9?  — 1 for  all  temperatures  [17,  34,  36,  37]. 

For  the  isothermal  case  with  T = Tm,  the  system  admits  a steady  one- dimensional 
solution  ip{x)  given  by 

if  r n?  11 

(16) 


^tanh 

X 

1 

[2\/2(^V^)J 

where  we  have  chosen  the  origin  so  that  (/^(O)  = 1/2.  This  solution  shows  that  the  width  of 
the  interfacial  layer  is  proportional  to  the  product  The  surface  free  energy  (4)  for  this 

solution  is  given  by 
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= T.e£ 


°°  2 . TM{(y/a) 


ip^dx  = 


6y/2a 


(17) 


Thus  to  maintain  a finite  surface  free  energy  in  the  sharp-interface  hmit  that  ( ^/a  tends  to 
zero  requires  the  constant  a to  tend  to  zero  as  well;  i.e.,  the  barrier  height  of  the  double- well 
potential  becomes  large. 


3 Anisotropic  Formulation 

To  describe  anisotropic  surface  free  energies,  we  allow  the  coefficient  ^ that  appears  in  the 
gradient  energy  term  ^^|Vy?|^/2  to  depend  on  the  orientation  0 of  the  contours  of  constant 
phase;  i.e.,  we  set  ( = ^(0)?  where 

0 = arctan((,cjy/(^a:)  (18) 

is  the  angle  that  the  normal  to  these  curves  makes  with  the  x-axis.  The  angle  0 is  thus  defined 
throughout  the  domain;  in  the  sharp-interface  limit  in  which  the  crystal-melt  interface  is 
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associated  with  the  curve  (p  = 1/2,  it  reduces  to  the  angle  9 that  the  interface  normal  makes 
with  a reference  direction.  The  appropriate  form  for  {(0)  that  is  necessary  to  recover  a given 
anisotropic  surface  free  energy  7(^)  will  be  noted  shortly.  Practical  difficulties  in  defining 
0 that  are  associated  with  the  fact  that  |V(^|  tends  to  zero  in  the  bulk  regions  far  from  the 
interface  are  inconsequential  to  the  asymptotic  analysis,  since  the  role  of  surface  free  energy 
is  significant  only  in  the  neighborhood  of  the  interface  where  |V(/?|  is  non-negUgible. 

The  term  in  Eq.  (8),  which  arises  from  the  variation 


if  ^ is  constant,  is  then  replaced  by  a more  complicated  term 


(19) 


~ + additional  terms. 


8(p  y 

whose  specific  form  we  next  compute. 
The  variation  of  the  integral 


(20) 


is  given  by  (here  d^jdQ) 

«/  = / dV 

= J {-^VV  • + ii'lfj'fiy  - dV 

= - / VV  + 2«'V0  • + (Kf  + a")  [<^X0„  - ¥>v0x]}  6v> dV, 


(21) 


(22) 


where  we  have  integrated  by  parts  and  discarded  the  boundary  term  because  it  does  not 
contribute  to  the  functional  derivative.  In  doing  this,  we  have  used  the  relation 

^x8<Py  - <PyS(p^ 


By  using  the  expression 


we  obtain 


6Q  = 


S/(p 


|V(^|  = 


cos  0 X + sin  0 y. 


|V¥>I 

Vt&v  - <fiy@x  = |Vv5|  V 


vivv’iy 


(23) 


(24) 


(25) 
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and 


V0.V^  = |VHz.Vx(^)i 


(26) 


here  x,  y,  and  z are  units  vectors  in  the  x,  y,  and  z directions,  respectively.  The  expression 
(20)  therefore  has  the  form 

- ^ (5[«0)f IVv^l^)  = K(0)pvV  + 2«0)a0)|Vv|z  Vx 

+ {\m)r + «0)r'(0))  iv^i  v ■ . (2?) 

An  additional  source  of  anisotropy  is  introduced  by  letting  the  empirical  relaxation  coeffi- 
cient depend  on  orientation  as  well,  so  that  t = t(0).  The  anisotropic  form  of  the  phase-field 
equations  (8)  and  (10)  then  becomes 


c{T,v)^  + L{T)ip'{v)^  = V • [WT] , 

■(0)g  = QiT)p'iv)  - G'M  - ± (iK(0)l^|Vvp| 


(28) 


(29) 


In  the  isothermal  case  with  T = Tm,  the  equations  admit  steady,  one- dimensional  solu- 
tions of  the  form  ip  = ip{x  • fi),  where  n is  a constant  unit  vector  and  x - n = s cos  6o  + ysin  Oq. 
The  orientation  is  then  constant,  with  0 = 6o,  and  the  solution  is  given  by  [c.f.  Eq.  (16)] 


■"=2 


tanh 


(. 


X • n 


+ 1 


\2V2C{eo)y^ 

The  surface  free  energy  for  this  orientation  is  then  given  by 

TM[({Go)\/a] 


(30) 


7(«o)  = 


6\/2a 


(31) 


This  one- dimensional  solution  shows  that  interface  width  also  varies  with  orientation;  this 
width  can  be  characterized  by  the  parameter 


p{eo)  = 6%/2[^(eo)\/a] 


(32) 


which  represents  the  width  of  the  transition  layer  from  (p  « 0.05  to  ~ 0.95.  We  note  that 
7(^0)  and  t]{9o)  are  both  proportional  to  ^(^o),  i-e.  they  have  the  same  anisotropy. 
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3.1  Dimensionless  Equations 


It  is  convenient  to  work  in  dimensionless  units.  We  choose  a length  scale  I that  represents 
a geometrical  length  in  the  system,  such  as  a typical  radius  of  curvature  of  the  macroscopic 
phase  boundaries.  We  choose  a diffusive  time  scale  C£,£^/A:£,,  where  ki,  is  the  hquid  thermal 
conductivity  and  cl  = deijdT  is  the  heat  capacity  per  unit  volume  of  the  liquid,  both 
evaluated  at  the  melting  point.  We  measure  temperature  relative  to  the  melting  point  in 
units  of  Tat,  and  measure  energy  density  in  units  of  L[Tm)- 

Appropriate  scalings  for  the  width  of  the  interfacial  layer  and  the  double-well  barrier 
height  are  incorporated  by  introducing  the  small  dimensionless  parameter  e,  defined  by 

o-L{Tm) 


e = 


Tm 


(33) 


where  L[Tm)  is  the  latent  heat  per  unit  volume  at  the  melting  point.  A thin  interfacial  layer 
is  obtained  by  setting 


= er(0). 


(34) 


where  r(0)  is  of  order  unity  and  is  a dimensionless  form  of  ^(0).  The  dimensionless  governing 
equations  may  then  be  written  in  the  form 


_du  1 - dip 

^-37+  oip(v>)^  = ^ 


dt 


dt 


kSJu 


(35) 


e^f 


(0)^  = -¥’)(</^’-1/2) + e(5(u)y((,o)-e2A  ^i[r(0)]2|Vyj|2^  , (36) 


where  the  space  and  time  variables  are  now  dimensionless,  and  the  variational  deriva- 
tive in  the  latter  expression  is  given  by  Eq.  (27)  with  r(0)  replacing  ^(0).  Here  u = 
[T  — Tm)/Tm  is  the  dimensionless  temperature,  k{u,ip)  = k{T,(p)lkL,  c(u,(p)  = c{T,(p)lcL, 
Liu)  = L{T)IL{Tm),  and 


Q{u)  = j 
Jo 


Tmcl 

L(TmY 

” L(u) 


S = 

f(0)  = 


(l+u)2 


du  = u -\-  O(u^). 


(37) 

(38) 

(39) 


9 


To  retain  the  effects  of  interface  kinetics,  f(0)  should  be  assumed  to  be  of  order  unity.  The 
dimensionless  Helmholtz  free  energy  functional  for  T = Tm  has  the  form 

From  Eqs.  (31)-(34),  we  see  that  the  dimensionless  surface  free  energy  for  orientation  0 = 

7(^0)  m) 

L{Tm)^  6v^ 

4 Matched  Asymptotic  Expansions 

An  asymptotic  expansion  of  the  dimensionless  phase-field  equations  (35)-(36)  in  the  sharp 
interface  Hmit  allows  the  identification  of  appropriate  forms  for  the  coefficients  ^(0),  t(0), 
and  a in  order  to  recover  the  anisotropic  generalization  (2)  of  the  Gibbs-Thomson  equation. 
The  formal  procedure  is  similar  to  that  employed  by  Caginalp  [29]  for  the  isotropic  case,  so 
we  provide  an  abbreviated  version  of  the  asymptotic  expansion.  To  perform  the  expansion, 
two  subregions  of  the  domain  are  identified:  in  the  inner  region,  which  represents  the  vicinity 
of  the  interfacial  layer,  the  gradient  of  the  phase  field  is  large  and  the  temperature  varies 
slowly,  and  in  the  outer  region,  which  represents  the  bulk  phases,  the  phase  field  is  essentially 
constant.  In  each  phase,  the  solution  can  be  represented  by  an  asymptotic  expansion  in  terms 
of  appropriately-scaled  inner  or  outer  variables.  The  inner  and  outer  regions  share  a common 
region  of  overlap,  and  in  this  intermediate  region,  the  asymptotic  expansions  for  the  inner 
and  outer  solutions  can  be  matched  to  determined  the  solution.  Roughly  speaking,  the  outer 
solution  determines  the  far-held  boundary  conditions  for  the  inner  solution,  and  the  inner 
solution  determines  the  appropriate  interfacial  jump  conditions  for  the  outer  solution. 

4.1  Outer  Solution 

The  outer  solution  is  dehned  in  the  bulk  phases  where  the  spatial  variation  of  is  small, 
and  this  variation  is  on  an  0(1)  length  scale.  The  solution  is  formally  expanded  in  powers 
of  e, 

u = -f  -f  e . . . , (42) 
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(1)  -L  .2,. (2) 


(43) 


The  leading  order  solution  is  given  by  = 1 and  = 0 in  the  liquid  and  solid  regions, 
respectively,  and  the  leading  order  temperature  in  each  region  satisfies  the  usual  diffusion 
equation 


:(o) 


~Tr 


= V • 


(0) 


(44) 


where  and  A:f°l  have  values  appropriate  to  the  respective  bulk  phases. 

Far-held  boundary  conditions  for  the  thermal  held  are  assumed  to  be  known,  but  interfa- 
cial jump  conditions,  that  hold  in  the  limit  that  the  interfacial  layer  becomes  sharp,  must  be 
determined  by  matching  with  the  inner  solution.  The  higher-order  corrections  for  thermal 
held  may  be  computed  by  continuing  the  procedure,  but  are  not  required  for  the  subsequent 
analysis.  From  Eq.  (36)  we  see  that  the  hrst-order  correction  for  the  phase  held  vanishes 
identically  under  the  assumption  that  p^(0)  = P^(l)  ==  0.  The  formal  expansion  for  the  outer 
solution  breaks  down  near  the  interfacial  layer,  where  the  variation  of  (/?  is  large. 


4.2  Inner  Solution 

To  perform  the  inner  expansion  in  the  interfacial  layer,  it  is  convenient  to  introduce  a local 
coordinate  system  based  on  a parametrization  of  the  curve  (p[xjy,t)  = 1/2.  In  terms  of 
the  arclength  s,  this  curve  may  be  expressed  in  the  form  x = X{s,t)  and  y = Y[s,t). 
The  curve  has  a tangent  vector  {X',Y'),  a normal  vector  {Y' ^ —X'),  and  a normal  velocity 
Vn  = Y'Xt  — X'Yt,  where  the  prime  denotes  the  derivative  with  respect  to  arclength  and 
time  derivatives  are  indicated  by  subscripts.  We  use  s as  one  of  the  local  coordinates,  and 
use  the  distance  r along  the  normal  as  the  other  coordinate,  so  that 

x{r,s,t)  = X{s,t) rY\s,t),  (45) 

y{r,s,t)  = Y{s,t)  -rX'{s,t).  (46) 

The  orientation  of  the  curve  is  chosen  so  that  the  soHd  lies  on  the  left  if  the  curve  is  traversed 
in  the  direction  of  increasing  s,  and  the  normal  then  points  into  the  liquid;  the  coordinate 
system  is  described  in  more  detail  in  the  Appendix.  The  governing  equations  in  the  local 
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coordinates  are  further  transformed  by  introducing  the  scaled  variable  p = rje  and  writ- 
ing ip{p,s,t)  = (p{x,y,t)  and  u{p,s,t)  = u{x,y,t)  in  the  inner  region,  with  corresponding 
expansions 


u = + • • • , 

(p  = -f  4-  + . . . 


(47) 

(48) 


in  terms  of  the  inner  variables. 

4.2.1  Matching  Conditions 

Matching  conditions  provide  the  far-held  boundary  conditions  for  the  inner  solution.  The 
outer  solution  u[x,y,t)  at  a point  near  the  curve  (/?  = 1/2  is  written  as  a function  of  the 
inner  variables,  and  the  resulting  expressions  are  expanded  in  e to  obtain 

n(X  + epY',Y-  epX\t)  = u^i\x,Y,t)  + £ |4''(X,r,4)  + p^(X,Y,t)\  + 0(^),  (49) 

where  the  plus  or  minus  subscript  indicates  the  limiting  behavior  of  the  outer  solution  as  p 
tends  to  zero  through  positive  or  negative  values,  respectively,  to  allow  for  the  possibility 
of  discontinuous  behavior  of  the  outer  solution  in  the  sharp- interface  limit.  The  limiting 
behavior  of  the  phase  held  is  simply  y?  = 1 -f-  O(e^)  for  p > 0 and  (,0  = (9(e^)  for  p < 0. 


4.2.2  Transformation  to  Inner  Variables 


The  orientation  angle  0,  which  in  the  inner  region  satishes 


, nr  ^ -{X'^epY")^,^eY'Cp, 

tan0(p,s,y?p,(^,)  = —— 


{Y'  - epX")(pp  -\-  eX'ifg 


can  be  expanded  in  the  form 


with 


0 = 0(0)  + e0(i)  + o{^). 


tan0(°)  = -X'jY'  - tan0(s). 


(50) 


(51) 


(52) 
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i.e.,  to  leading  order,  0^°)  is  simply  the  normal  angle  0(s)  to  the  interface,  and  is  independent 
of  the  variable  p.  A short  calculation  shows  that  the  first  order  correction  is  given  by 


0(1)  = ^(0)/^(0). 


(53) 


we  note  that  s)  vanishes  when  p = 0,  but  is  not  necessarily  zero  if  p ^ 0. 

The  governing  equations  may  be  transformed  by  using  the  results  given  in  the  Appendix; 
here  we  note  some  of  the  intermediate  expressions  before  presenting  the  final  results.  The 
time  derivative  transforms  according  to 

dip  —Vri  dip^^^ 


dt  e dp 

The  transformed  Laplacian  assumes  the  form 

1 


+ 0(1). 


(54) 


vv  = 


(1  + cpK) 


[(1  + C/)A')y?r]r  + 
K(s) 


(1  + epK) 


= i(^W  + e,?W)  + :^^(0)  + o(i), 
and  the  gradient  is  given  by 

|V<^|  = + + 0(e); 


(55) 


(56) 


here  we  assume  that  the  coordinates  are  oriented  so  that  > 0 when  simphfying  the 
square  root. 

We  also  have  the  expansions 


1 

— < 

1 

(1  + epK)<pp 

+ e 

p 

ipj{l  + epK) 

> 

> 

(1  + epK) 

e 

+ eV^/(l  + epX)2_ 

3 > 

- K{s)  + 0{e), 


(57) 


and 


z • V X 


/V^\  _ 1 1 

' 

1 

^p 

> 

(l  + e/pif) 

e 

+ eVa/Cl  + 

p 

+ eV^/(l  + epir)2_ 

a > 

+ 0(e). 


(58) 
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We  therefore  have 


ir(e)l^<?JS>  + € {[r(0)]^<?w  + [r(0)i^A'(5)v('’)  + 2r(0)r'(0) 


X 


(3(0)  _ ^PP 
(3(0) 
<Pp 


+ ([r(0)f  + r(0)r"(0))  + o(e^) 


+ ([r'(e)]2  + r(e)r"(e))  } + o{e%  (59) 


where  in  the  final  expression  we  have  expanded  F = r(0)  using 

0 9{s)  + + O(e^). 

to  simplify  the  result. 


(60) 


4.2.3  Leading-Order  Solution 


The  leading-order  equations  take  the  form 


[r(«)]^^  - = 0, 


dp^ 


where 


1 


5(^(0))  = - 1)(<^»)  - 1/2). 


(61) 

(62) 

(63) 


Here  denotes  the  thermal  conductivity  evaluated  with  the  leading  order 

solution. 

The  thermal  matching  conditions  for  the  inner  problem, 

u>-'>\p,s,t)-*u<i\X,Y,t)  (64) 


as  p — > ±00  implies  that  j d p = 0,  and  so  is  independent  of  p.  It  follows  that 

the  leading  order  thermal  field  for  the  outer  solution  is  continuous  across  the  interface. 

The  leading  order  solution  for  the  phase  field  is  [cf.  Eq.  (16)] 

^(%,.)  = i{tanh(j^)+l);  (65) 

here  the  5-dependence  of  enters  through  the  function  F,  whose  argument  is  the  normal 
angle  9(s)  along  the  interface. 
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4.2.4  First-Order  Solution 


The  first-order  inner  problem  has  the  form 


(66) 


- 2r(e)r'(ff).?<,°>  - (|r'(«)]^  + r(«)r"(fl))  (67) 


The  thermal  matching  conditions  for  the  inner  problem  imply  that 

as  p — > ±oo.  Integration  of  the  thermal  equation  then  gives  the  equation 


(68) 


- V, 


I — — = f^L—^ «5‘ 


dn 


dn 


(69) 


which  is  the  appropriate  heat  flux  boundary  condition  for  the  outer  solution. 

Differentiation  of  Eq.  (62)  with  respect  to  p shows  that  the  function  is  a homogeneous 
solution  of  Eq.  (67).  The  right  hand  side  of  this  equation  must  be  then  be  orthogonal  to  this 
function,  which  provides  the  solvability  condition 

/OO  /•OO 

(<?(»))^  dp  = -Qiu^o))  - [r(«)]^A'(s)  / dp 

-OO  j — OO 

f]  r roo  *]  / . TOO 

- r(fi)r'(9)-  dp\  - ([r'(e)]^  + r(9)r"(«))  k(s)  dp,  (70) 

or,  using  the  relation 

’)  dp  - gy2r(e)’ 


(71) 


we  have 


^ = -W<»>)-^ir(-i)  + rw]x(.), 


(72) 


which  is  our  principal  result. 
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5 Discussion 


The  expression  (72)  is  the  dimensionless  version  of  Eq.  (2)  if  we  take  Q{u)  « u and  revert 
to  dimensional  variables, 

Tm 


^{6) 


= Tm  — 


L{Tm) 


(7  + lee)K  — Ti\ 


here  we  have  identified 


^l(e) 


_ (&s/2L(Tm)\  m-Ji 


\ 


rp2 

-‘-M 


j r{») 


(73) 


(74) 


as  the  appropriate  form  for  the  kinetic  coefficient,  and  we  recall  that  the  surface  free  energy 
is  given  by 


7(0)  - 


6\/2a 


(75) 


Note  that  anisotropy  in  the  parameter  ((6)  induces  anisotropy  in  both  the  surface  free 
energy  and  the  kinetic  coefficient,  even  when  the  parameter  r is  isotropic;  note  also  that  the 
anisotropy  of  the  kinetic  coefficient  depends  on  the  ratio  of  ( and  r. 

The  expressions  (74),  (75),  and  (32)  relate  the  physical  parameters  /x,  7,  and  the  interface 
width  7)  to  the  phase-field  parameters  r,  and  a.  These  relations  can  be  inverted  to  yield 
expressions  for  the  phase- field  parameters  in  terms  of  the  physical  parameters,  viz. 


L(T„),(0) 

W)  ’ 

(76) 

J-M 

(77) 

Tm7(0) 

""  " 727(0)  ’ 

(78) 

note  that  the  latter  expression  impfies  a common  functional  dependence  for  rj  and  7 in  order 
for  the  resulting  parameter  a to  be  constant.  Alternatively,  the  parameters  r and  ^ may  be 
expressed  in  terms  of  7,  fi,  and  the  single  parameter  a by  using  Eq.  (78)  to  eliminate  rj  in 
equations  (76)  and  (77). 

As  mentioned  in  the  introduction,  an  anisotropic  surface  free  energy  can  also  be  obtained 
by  replacing  the  isotropic  gradient  energy  term  by  the  more  general  quadratic  form  -|- 

^yiPy)/2  [23],  where  and  (y  are  constants.  In  this  case,  the  term  in  Eq.  (8)  would 
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be  replaced  by 


- ^ ff  + I (79) 

If  X is  scaled  by  (x  a-nd  V is  scaled  by  (y,  this  expression  reduces  to  that  for  the  isotropic  case. 
Since  the  isotropic  equihbrium  shape  is  given  by  a circle,  the  corresponding  anisotropic  equi- 
librium shape  is  therefore  an  elhpse.  This  choice  of  gradient  energy  leads  to  one-dimensional 
solutions  of  the  form  (30)  with  surface  free  energy  (31),  where 

^{do)  = yjil  C0S2  do  + ^2  siji2 

This  particular  form  of  gradient  energy  term  thus  leads  to  a surface  free  energy  with  two-fold 
axes  of  symmetry  about  the  orientations  0 = 0 and  6 = 7r/2,  and  possessing  only  two  degrees 
of  freedom  (^^  and  ^y).  The  present  approach  for  introducing  anisotropy,  wherein  a scalar 
function  ( = ^(0)  is  instead  employed,  allows  general  surface  energies  7 = 7(0)  to  be  treated. 
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Appendix 

We  collect  some  pertinent  facts  about  the  local  coordinate  system  used  in  the  inner  expansion 
of  the  matching  procedure.  The  orthogonal  coordinates  r and  s are  defined  relative  to  the 
moving  curve  x = Ar(5,i)  and  y = Y(s,t),  where  s is  arclength  along  the  curve.  The 


coordinate  transformation  is  given  by 

x(r,5,t)  = 

X{s,t)^rY'{s,t) 

(81) 

y{r,s,t)  = 

Y[s,t)  — rX'{s,t). 

(82) 
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Thus 


X. 


(83) 

(84) 


Xr  = Y’{s,t), 
X,  = X'(s,t)  + rY"{s,t), 


Vt  = -X'{s,t), 

y.  = Y'{s,t)~rX"(s,t). 


If  the  angle  that  the  normal  to  the  interface  makes  with  respect  to  the  x-axis  is  denoted  by 
^(5),  then  6g  = /C(s)  is  the  local  curvature  of  the  interface.  We  have 

X‘  + iV  = (85) 


and  by  differentiating  we  have 


X"  + iY"  = -e^^K. 


It  follows  that 


and 


K = -{X"  - iY")e^^  = X'Y"  - Y'X", 


{X"f  + {Y''f  = K\ 


For  the  unit  circle,  this  expression  gives  a positive  curvature  K = 1. 
We  have  the  Jacobian 


(86) 

(87) 

(88) 


h{r,s)  = XrV,  - x.yr  = 1 +rK{s).  (89) 

Since  XrXg  + = 0,  the  coordinates  are  orthogonal,  and  the  square  of  the  element  of 

differential  arclength  dS  in  the  three-dimensional  set  of  orthogonal  coordinates  {r,s,z)  is 
given  by 

dS^  = dr^  h?ds^  -|-  dz^ . (90) 

The  gradient  of  a function  transforms  according  to 

VV?  = V’r?  + 7V’5S,  (91) 

h 

and  the  Laplacian  is  given  by 

VV  = i{(AW.  + (iV'.)j.  (92) 

Given  a vector 

A(r,  5)  ==  u(r,  5)? -f  t;(r,  5)s,  (93) 
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we  have 


(94) 


V ■ A = i{(/iu)r  + t),}, 

tl 

and 

z • V X A = T{i^'^)r  — (95) 

h 

(These  relations  are  analogous  to  those  in  cylindrical  coordinates,  if  h[r,  s)  is  identified  with 
the  radius  in  the  corresponding  expressions.) 

If  we  regard  the  coordinates  r and  s as  functions  of  x,  y,  and  t,  we  note  that  inverting 
the  Jacobian  matrix  gives  the  relations 


hr^  = Y'{s)  - rX"{s), 

hr,  = -X‘{s)-rY''{s), 

(96) 

hs,  = X'{s), 

h,y  = r(s). 

(97) 

We  also  have 

hrt  = -{XtY'  - YtX')  + r{XtX"  + YY") 

+ r^{Y;X"  - X^Y")  = -Vn  + 0(r),  (98) 

hst  = -{XtX'  + YtY')  + t{Y'X[  - XX) 

= — ntan  + 0{r),  (99) 

where  Vn  = Y'Xt  — X'Yt  is  the  normal  velocity  of  the  curve  and  Utan  = Y'Yt  + X'Xt  is  a 
tangential  velocity  which  depends  on  the  specific  choice  of  arclength  parametrization. 
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